*
* $Id$
*
* $Log: gnopg1.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:39  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:17  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:52  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
*-- Author :
      SUBROUTINE GNOPG1(X,P,SNXT)
*********************************************************************
***** GNOPG1 ********************************************************
*
*     GNOPG1 ... 20-DEC-1988
*     Version 1.1
*     Rolf Nierhaus
*
*********************************************************************
*
*     Copyright   CERN,   Geneva  1988  -  Copyright  and  any  other
*     appropriate legal protection of  these  computer  programs  and
*     associated  documentation  reserved  in  all  countries  of the
*     world.
*
*********************************************************************
*
*          Subroutine  GNOPG1 is called by GNOPGO for the computation
*     of SNXT, the distance from a point P  along  a  track  T  to  a
*     boundary  surface  of a Geant volume V of shape PGON. The point
*     P is outside the volume V. If the track T has  no  intersection
*     with the volume V, the vector distance SNXT is set to 1.E10.
*
*          V  is  generally  a composite volume consisting of several
*     sections. The sections have  boundary  surfaces  orthogonal  to
*     the   Z-axis.   Each  section  consists  generally  of  several
*     sectors. Each sector is an  "elementary"  convex  volume.  This
*     package  assumes it is either a hexahedron or a pentahedron. If
*     it is a pentahedron, it has 6 vertices, of  which  two  are  on
*     the  Z-axis.  All  sectors  of  the same section are congruent.
*     Each section has the same number of sectors.
*
*          GNOPG1 calls GNOPG2 for each section, GNOPG2  call  GNOPG3
*     for  each  sector.  GNOPG4  is  called  to  store  the  surface
*     parameters of  a  sector  in  the  common  block  GCQ1.  GNOPG6
*     computes  the  vector  distance  to  a  convex  volume.  GNOPG7
*     computes the  vector  distance  to  a  plane  surface.  Logical
*     function  GNOPG8  determines  if  a  point  is  inside a convex
*     volume, and logical function GNOPG9 determines if  a  point  is
*     inside a region delimited by a plane surface.
*
*          We  describe each surface by 6 parameters: the first three
*     are   the   coordinates   of   a   point   on    the    surface
*     XS(I),YS(I),ZS(I),  the  other  three are the components of the
*     normal vector of the surface XN(I),YN(I),ZN(I). I is the  index
*     of  the surface. We consider only one sector at a time, and the
*     number of boundary  surfaces  is  never  larger  then  6.  Each
*     surface  divides  the  space  into  two  regions:  the positive
*     region and the negative region. We choose the direction of  the
*     normal  vectors  of the boundary surfaces such that the bounded
*     volume is within the positive region of each surface, that  is,
*     the normal vector is pointing to the inside of the volume.
*
*          Logical   function   GNOPG9  returns  TRUE  if  the  point
*     (XP,YP,ZP) is within the positive region of  the  surface  with
*     index   I.   This   is  the  case  if  the  scalar  product  of
*     (XP-XS,YP-YS,ZP-ZS) and (XN,YN,ZN) is positive (or zero).
*
*          GNOPG8 is true if  the  point  (XP,YP,ZP)  is  within  the
*     positive region of all bounding surfaces.
*
*          To  find  the  distance  from  a  point (XP,YP,ZP) along a
*     track  with  direction  cosines   (XD,YD,ZD)   to   a   surface
*     (XS,YS,ZS)(XN,YN,ZN),  we  compute  first the scalar product of
*     the  vector  (XS-XP,YS-YP,ZS-ZP)   with   the   normal   vector
*     (XN,YN,ZN),  then  the scalar product of the vectors (XD,YD,ZD)
*     and (XN,YN,ZN).  The  first  scalar  product  is  the  shortest
*     distance  from  the  point  to  the  plane,  the  second scalar
*     product is the cosine of the angle between the  track  and  the
*     plane  normal.  The  quotient  is  the vector distance. If this
*     vector distance is  positive  (or  zero)  we  set  the  logical
*     variable  FLAG  TRUE.  GNOPG7  is  called with three parameters
*     I,FLAG and DIST. I is the index of the  surface,  and  DIST  is
*     the vector distance if FLAG is TRUE.
*
*          To  find the vector distance from a point to an elementary
*     volume, all bounding surfaces of the volume are considered.  If
*     the  point  is  in  the positive region of a surface, the track
*     could possibly exit through the surface, but  it  cannot  enter
*     through  it.  For  those surfaces which have the point in their
*     negative region, we  determine  if  the  track  intersects  the
*     surface,  and  what  is the distance to the intersection point.
*     Only the largest of these distances  is  a  candidate  for  the
*     distance  from  the point to the volume. It is however possible
*     that the track misses the elementary volume entirely.  This  we
*     find  out  by applying the function GNOPG8 (Inside volume test)
*     to the coordinates of the intersection point. GNOPG6 is  called
*     with  two parameters: a logical variable FLAG, which is TRUE if
*     the track  intersects  the  volume,  and  DIST,  which  is  the
*     distance if FLAG is TRUE.
*
*          The  coordinates of  the point P and the direction cosines
*     of the  track  T  are  stored  in  the  common  block  GCQ2  by
*     subroutine  GNOPG1. The parameter array P in the call to GNOPG1
*     contains in its first two members the lower  phi-limit  of  the
*     PGON  and  the range in phi. Both angles are in degrees. GNOPG1
*     convertes from degrees to radians and stores the phi-limits  of
*     the  first  sector  of  each  section  in the common block GCQ3
*     together with the number of sectors per section. The number  of
*     sectors  per  section  is  contained in the third member of the
*     parameter array P. The other members of P  have  the  following
*     meaning:  From  P(5)  onwards,  there are groups of three. Each
*     group describes a section boundary. The  first  member  is  the
*     Z-coordinate  of the boundary. The second and the third are the
*     distances from the Z-axis to the inner and outer edges at  that
*     boundary.  If  there is no inner edge, the sector is limited by
*     the Z-axis. In this case the second members of  the  group  are
*     zero.  Groups may be shared by adjacent sections. Otherwise the
*     Z-coordinates of two consecutive  groups  are  the  same.  P(4)
*     contains the number of groups.
*
*          GNOPG1  calls  GNOPG2  with  8 parameters. The first 6 are
*     input parameters, the first  3  refer  to  the  "left"  section
*     boundary,  the next 3 to the "right" section boundary. The last
*     two parameters are output. Logical variable  FLAG  is  TRUE  if
*     the  track  intersects  the  section.  In this case DIST is the
*     distance.
*
*          GNOPG2 calls GNOPG3 with four parameters. The first 2  are
*     input  parameters,  namely  the  phi-limits  of the sector. The
*     last two parameters again are  output:  FLAG  is  TRUE  if  the
*     track   intersects  the  sector.  In  this  case  DIST  is  the
*     distance. The other variables required by GNOPG3 are  the  same
*     for  all  sectors  of the same section and are stored by GNOPG2
*     in the common block GCQ4. GNOPG2  also  stores  in  the  common
*     block GCQ5 the number of boundary surfaces of each sector.
*
*          The  surfaces  orthogonal  to  the Z-axis are the same for
*     all sectors of a section. The surface parameters of  these  two
*     sections  are  stored  by  calls  to  GNOPG4  from  GNOPG2. The
*     surface parameters of the other  three  or  four  surfaces  are
*     stored from GNOPG3.
*
*          GNOPG3  sets  FLAG  TRUE,  if  the  track T intersects the
*     corresponding sector. GNOPG2 finds  the  shortest  distance  to
*     all  intersected  sectors,  and  set  FLAG TRUE, if the track T
*     intersects  the  corresponding  section.   GNOPG1   finds   the
*     shortest  distance  to  all intersected sections. If no section
*     intersects, the track FLAG is set FALSE, and 1.E10 is  returned
*     as SNXT.
*
***** Subroutine GNOPG1 *************************** 04-DEC-1988 *****
#if !defined(CERNLIB_SINGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (F=0.01745329251994330D0)
#endif
#if defined(CERNLIB_SINGLE)
      PARAMETER (F=0.01745329251994330)
#endif
      REAL X(6),P(49),SNXT
      LOGICAL FLAG,FLAG1
      DIMENSION XS(6), YS(6), ZS(6), XN(6), YN(6), ZN(6)
      LOGICAL FLAG1X, GNOP1X, GNOP2X, GNOP4X
      PARAMETER (ONE=1,HALF=ONE/2)
      XP=X(1)
      YP=X(2)
      ZP=X(3)
      XD=X(4)
      YD=X(5)
      ZD=X(6)
      NT=P(3)+.5
      P1=F*P(1)
      P2=F*P(2)/NT
      INDEX=5
      MINDEX=3.*P(4)+1.5
      FLAG=.FALSE.
      DIST=1.E10
   10 Z1=P(INDEX)
      D1N=P(INDEX+1)
      D1X=P(INDEX+2)
      Z2=P(INDEX+3)
      D2N=P(INDEX+4)
      D2X=P(INDEX+5)
C*****  Code Expanded From Routine:  GNOPG2
C*****  Code Expanded From Routine:  GNOPG4
      XS(1) = 0.
      YS(1) = 0.
      ZS(1) = Z1
      XN(1) = 0.
      YN(1) = 0.
      ZN(1) = 1.
C*****  End of Code Expanded From Routine:  GNOPG4
C*****  Code Expanded From Routine:  GNOPG4
      XS(2) = 0.
      YS(2) = 0.
      ZS(2) = Z2
      XN(2) = 0.
      YN(2) = 0.
      ZN(2) = -1.
C*****  End of Code Expanded From Routine:  GNOPG4
      P3 = P1
      P4 = P1 + P2
      Z3 = HALF*(Z1 + Z2)
      D3X = HALF*(D1X + D2X)
      TH1 = ATAN((D2X - D1X)/(Z2 - Z1))
      COSTH1 = COS(TH1)
      SINTH1 = SIN(TH1)
      D3N = HALF*(D1N + D2N)
      ISMAX = 5
      IF (D3N .NE. 0.) THEN
         ISMAX = 6
         TH2 = ATAN((D2N - D1N)/(Z2 - Z1))
         COSTH2 = COS(TH2)
         SINTH2 = SIN(TH2)
      ENDIF
      FLAG1 = .FALSE.
      DIST1 = 1.E10
      DO 60  I = 1, NT
C*****  Code Expanded From Routine:  GNOPG3
C*****  Code Expanded From Routine:  GNOPG4
         XS(3) = 0.
         YS(3) = 0.
         ZS(3) = Z3
         XN(3) = -SIN(P3)
         YN(3) = COS(P3)
         ZN(3) = 0.
C*****  End of Code Expanded From Routine:  GNOPG4
C*****  Code Expanded From Routine:  GNOPG4
         XS(4) = 0.
         YS(4) = 0.
         ZS(4) = Z3
         XN(4) = SIN(P4)
         YN(4) = -COS(P4)
         ZN(4) = 0.
C*****  End of Code Expanded From Routine:  GNOPG4
         P1X = HALF*(P3 + P4)
         COSP = COS(P1X)
         SINP = SIN(P1X)
C*****  Code Expanded From Routine:  GNOPG4
         XS(5) = D3X*COSP
         YS(5) = D3X*SINP
         ZS(5) = Z3
         XN(5) = -COSP*COSTH1
         YN(5) = -SINP*COSTH1
         ZN(5) = SINTH1
C*****  End of Code Expanded From Routine:  GNOPG4
         IF (D3N .NE. 0.) THEN
C*****  Code Expanded From Routine:  GNOPG4
            XS(6) = D3N*COSP
            YS(6) = D3N*SINP
            ZS(6) = Z3
            XN(6) = COSP*COSTH2
            YN(6) = SINP*COSTH2
            ZN(6) = -SINTH2
C*****  End of Code Expanded From Routine:  GNOPG4
         ENDIF
C*****  Code Expanded From Routine:  GNOPG6
         FLAG1X = .FALSE.
         DIST2X = 0.
         DO 20  IS = 1, ISMAX
C*****  Code Expanded From Routine:  GNOPG9
*     TRUE if (XP,YP,ZP) in positive region of surface I
            GNOP1X = 0. .LE. (XP-XS(IS))*XN(IS)+(YP-YS(IS))*YN(IS)+(ZP-
     +         ZS(IS))*ZN(IS)
C*****  End of Code Expanded From Routine:  GNOPG9
C*****  Code Expanded From Routine:  GNOPG7
            IF (.NOT.GNOP1X) THEN
               SPPMSN = (XP - XS(IS))*XN(IS) + (YP - YS(IS))*YN(IS) + (
     +            ZP - ZS(IS))*ZN(IS)
               SPDN = XD*XN(IS) + YD*YN(IS) + ZD*ZN(IS)
               IF (SPDN .EQ. 0.) THEN
                  DIST1X = -.000001
               ELSE
                  DIST1X = -SPPMSN/SPDN
               ENDIF
C*****  End of Code Expanded From Routine:  GNOPG7
               IF ((-1.E-5) .LE. DIST1X) THEN
                  FLAG1X = .TRUE.
                  DIST2X = MAX(DIST1X,DIST2X)
               ENDIF
            ENDIF
   20    CONTINUE
         IF (.NOT.FLAG1X) GO TO 50
         T = DIST2X + .001
         XQ = XP + T*XD
         YQ = YP + T*YD
         ZQ = ZP + T*ZD
C*****  Code Expanded From Routine:  GNOPG8
*     TRUE if (XP,YP,ZP) in volume
         GNOP2X = .FALSE.
         DO 30  IS1X = 1, ISMAX
C*****  Code Expanded From Routine:  GNOPG9
*     TRUE if (XP,YP,ZP) in positive region of surface I
            GNOP4X = 0. .LE. (XQ-XS(IS1X))*XN(IS1X)+(YQ-YS(IS1X))*YN(
     +         IS1X)+(ZQ-ZS(IS1X))*ZN(IS1X)
            IF (.NOT.GNOP4X) GO TO 40
C*****  End of Code Expanded From Routine:  GNOPG9
   30    CONTINUE
         GNOP2X = .TRUE.
   40    CONTINUE
         FLAG1X = GNOP2X
C*****  End of Code Expanded From Routine:  GNOPG8
   50    CONTINUE
C*****  End of Code Expanded From Routine:  GNOPG3
         IF (FLAG1X) THEN
            FLAG1 = .TRUE.
            DIST1 = MIN(DIST2X,DIST1)
         ENDIF
         P3 = P4
         P4 = P4 + P2
   60 CONTINUE
C*****  End of Code Expanded From Routine:  GNOPG2
      IF (FLAG1) THEN
         FLAG=.TRUE.
         IF (DIST1.LT.DIST) DIST=DIST1
      END IF
      INDEX=INDEX+3
      IF (MINDEX.LT.INDEX) THEN
         IF(FLAG) THEN
            SNXT=DIST
         ELSE
            SNXT=1.E10
         ENDIF
      ELSE
         IF (P(INDEX+3).EQ.Z2) INDEX=INDEX+3
         GO TO 10
      ENDIF
      END
